#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
睡眠时间预测系统

此模块实现了一个完整的睡眠时间预测系统，通过分析生活方式因素（如锻炼时间、阅读时间、手机使用时间等）
来预测个体的睡眠时间。系统包含数据预处理、特征工程、模型训练、超参数优化和结果分析等完整流程。

作者: AI助手
日期: 2024
版本: 1.0
"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_score
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.svm import SVR
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
import joblib
import warnings

# 忽略警告信息，保持输出简洁
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

class SleepPredictionSystem:
    """
    睡眠时间预测系统 - 用于分析和预测睡眠时间的机器学习模型
    
    该类实现了一个端到端的机器学习系统，用于基于各种生活方式因素预测睡眠时间。
    系统包含完整的数据处理流水线，从数据加载、预处理、特征工程到模型训练和评估。
    
    属性:
        df: 原始数据DataFrame
        X: 特征矩阵
        y: 目标变量（睡眠时间）
        scaler: 特征标准化器
        best_model: 性能最佳的模型
        best_model_name: 最佳模型名称
        results: 模型评估结果
        X_train, X_test, y_train, y_test: 训练集和测试集数据
    """
    
    def __init__(self):
        """
        初始化预测系统
        
        初始化系统所需的所有属性，包括数据存储、模型存储和评估结果等。
        """
        # 数据相关属性
        self.df = None  # 原始数据
        self.original_shape = None  # 原始数据形状
        self.X = None  # 特征矩阵
        self.y = None  # 目标变量（睡眠时间）
        
        # 数据集分割
        self.X_train = None  # 训练集特征
        self.X_test = None  # 测试集特征
        self.y_train = None  # 训练集目标变量
        self.y_test = None  # 测试集目标变量
        
        # 特征处理
        self.scaler = StandardScaler()  # 特征标准化器
        self.X_train_scaled = None  # 标准化后的训练集特征
        self.X_test_scaled = None  # 标准化后的测试集特征
        
        # 模型和评估
        self.best_model = None  # 性能最佳的模型
        self.best_model_name = None  # 最佳模型名称
        self.results = {}  # 模型评估结果字典
    
    def load_data(self, file_path='sleep.csv', encoding='gb2312'):
        """
        加载和预处理数据
        
        从CSV文件中读取睡眠时间预测数据集，重命名列名为英文，并显示数据集的基本信息。
        
        参数:
            file_path (str): 数据文件路径，默认为'sleep.csv'
            encoding (str): 文件编码，默认为'gb2312'以支持中文
            
        返回:
            self: 返回实例本身以支持方法链式调用
        """
        print("=== 加载和预处理数据 ===")
        
        # 加载数据
        self.df = pd.read_csv(file_path, encoding=encoding)
        self.original_shape = self.df.shape
        
        # 重命名列名为英文，便于后续处理和模型训练
        self.df.columns = ['exercise_time', 'reading_time', 'phone_time', 'work_time', 
                          'caffeine_intake', 'relax_time', 'sleep_time']
        
        print(f"原始数据集形状: {self.original_shape}")
        print("数据前5行:")
        print(self.df.head())
        
        return self
    
    def handle_outliers(self, method='iqr'):
        """
        处理异常值
        
        使用指定的方法检测并处理数据中的异常值，避免异常值对模型训练产生负面影响。
        
        参数:
            method (str): 异常值处理方法，可选值为'iqr'(四分位距法)或'zscore'(Z分数法)，默认为'iqr'
            
        返回:
            self: 返回实例本身以支持方法链式调用
        """
        print("\n=== 异常值处理 ===")
        
        # 创建数据副本以避免修改原始数据
        df_clean = self.df.copy()
        
        outliers_info = {}  # 存储每个特征的异常值数量
        
        for col in df_clean.columns:
            if method == 'iqr':
                # IQR方法 - 计算四分位数和四分位距
                Q1 = df_clean[col].quantile(0.25)  # 第一四分位数
                Q3 = df_clean[col].quantile(0.75)  # 第三四分位数
                IQR = Q3 - Q1  # 四分位距
                lower_bound = Q1 - 1.5 * IQR  # 下界
                upper_bound = Q3 + 1.5 * IQR  # 上界
                
                outliers = df_clean[(df_clean[col] < lower_bound) | (df_clean[col] > upper_bound)]
                outliers_info[col] = len(outliers)
                
                # 替换异常值为边界值，保持数据分布的合理性
                df_clean[col] = np.where(df_clean[col] < lower_bound, lower_bound, df_clean[col])
                df_clean[col] = np.where(df_clean[col] > upper_bound, upper_bound, df_clean[col])
            
            elif method == 'zscore':
                # Z-score方法 - 计算每个数据点的标准化得分
                z_scores = np.abs((df_clean[col] - df_clean[col].mean()) / df_clean[col].std())
                outliers = df_clean[z_scores > 3]  # 通常认为Z分数大于3的为异常值
                outliers_info[col] = len(outliers)
                
                # 替换异常值为均值
                df_clean[col] = np.where(z_scores > 3, df_clean[col].mean(), df_clean[col])
        
        # 更新数据
        self.df = df_clean
        
        print(f"异常值处理后数据集形状: {self.df.shape}")
        print("各特征异常值数量:")
        for col, count in outliers_info.items():
            print(f"  {col}: {count} 个")
        
        return self
    
    def exploratory_data_analysis(self):
        """
        探索性数据分析
        """
        print("\n=== 探索性数据分析 ===")
        
        # 基本统计
        print("\n数据描述性统计:")
        print(self.df.describe())
        
        # 缺失值检查
        print("\n缺失值检查:")
        missing_values = self.df.isnull().sum()
        if missing_values.sum() > 0:
            print(missing_values)
            # 如果有缺失值，填充
            self.df = self.df.fillna(self.df.mean())
            print("已填充缺失值")
        else:
            print("无缺失值")
        
        # 相关性分析
        print("\n=== 相关性分析 ===")
        correlation_matrix = self.df.corr()
        
        # 创建相关性热力图
        plt.figure(figsize=(12, 8))
        sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', center=0, 
                   square=True, linewidths=0.5, fmt='.2f')
        plt.title('生活方式因素与睡眠时间的相关性热力图', fontsize=16)
        plt.tight_layout()
        plt.savefig('相关性热力图.png', dpi=300, bbox_inches='tight')
        plt.show()
        
        # 打印与睡眠时间的相关性
        print("\n与睡眠时间的相关性（从高到低）:")
        sleep_correlations = correlation_matrix['sleep_time'].sort_values(ascending=False)
        for factor, corr in sleep_correlations.items():
            if factor != 'sleep_time':
                print(f"{factor}: {corr:.3f}")
        
        # 绘制特征与目标变量的散点图
        self._plot_feature_scatter()
        
        # 绘制目标变量分布
        self._plot_target_distribution()
        
        return self
    
    def _plot_feature_scatter(self):
        """绘制特征与目标变量的散点图"""
        fig, axes = plt.subplots(2, 3, figsize=(15, 10))
        axes = axes.ravel()
        
        features = ['exercise_time', 'reading_time', 'phone_time', 'work_time', 'caffeine_intake', 'relax_time']
        feature_names = ['运动时间', '阅读时间', '手机使用时间', '工作时间', '咖啡因摄入量', '放松时间']
        
        for i, (feature, name) in enumerate(zip(features, feature_names)):
            axes[i].scatter(self.df[feature], self.df['sleep_time'], alpha=0.6, s=20)
            axes[i].set_xlabel(name)
            axes[i].set_ylabel('睡眠时间')
            axes[i].set_title(f'{name} vs 睡眠时间')
            
            # 添加趋势线
            z = np.polyfit(self.df[feature], self.df['sleep_time'], 1)
            p = np.poly1d(z)
            axes[i].plot(self.df[feature], p(self.df[feature]), "r--", alpha=0.8)
        
        plt.tight_layout()
        plt.savefig('散点图矩阵.png', dpi=300, bbox_inches='tight')
        plt.show()
    
    def _plot_target_distribution(self):
        """绘制目标变量分布"""
        plt.figure(figsize=(14, 6))
        
        plt.subplot(1, 2, 1)
        plt.hist(self.df['sleep_time'], bins=30, alpha=0.7, color='skyblue', edgecolor='black')
        plt.xlabel('睡眠时间 (小时)')
        plt.ylabel('频次')
        plt.title('睡眠时间分布直方图')
        plt.grid(True, alpha=0.3)
        
        plt.subplot(1, 2, 2)
        plt.boxplot(self.df['sleep_time'])
        plt.ylabel('睡眠时间 (小时)')
        plt.title('睡眠时间箱线图')
        plt.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.savefig('睡眠时间分布.png', dpi=300, bbox_inches='tight')
        plt.show()
        
        print(f"\n睡眠时间统计信息:")
        print(f"平均值: {self.df['sleep_time'].mean():.2f} 小时")
        print(f"中位数: {self.df['sleep_time'].median():.2f} 小时")
        print(f"标准差: {self.df['sleep_time'].std():.2f} 小时")
        print(f"最小值: {self.df['sleep_time'].min():.2f} 小时")
        print(f"最大值: {self.df['sleep_time'].max():.2f} 小时")
    
    def feature_engineering(self):
        """
        特征工程 - 创建更多有意义的特征
        
        此方法通过组合原始特征创建新的衍生特征，增强模型的预测能力:
        1. 准备基础特征矩阵和目标变量
        2. 生成衍生特征来捕捉变量间的关系和复杂模式
        3. 分析新特征与目标变量的相关性
        
        衍生特征解释:
        - total_activities: 日常活动总量，包括运动、阅读、手机使用和放松时间
        - work_life_balance: 工作与生活平衡指标，衡量工作时间相对于放松活动的比例
        - caffeine_per_hour: 每工作小时的咖啡因摄入量
        - screen_time_ratio: 屏幕时间占总活动时间的比例
        - healthy_activity_index: 健康活动指数，衡量健康活动与屏幕时间的比率
        - daily_time_utilization: 每日时间利用总和
        """
        print("\n=== 特征工程 ===")
        
        # 准备特征和目标变量
        self.X = self.df.drop('sleep_time', axis=1)
        self.y = self.df['sleep_time']
        
        # 创建衍生特征
        # 1. 活动总量 - 所有日常活动的总时间
        self.X['total_activities'] = (self.X['exercise_time'] + self.X['reading_time'] + 
                                     self.X['phone_time'] + self.X['relax_time'])
        
        # 2. 工作生活平衡指数 - 较低的值表示更好的平衡
        self.X['work_life_balance'] = self.X['work_time'] / (self.X['exercise_time'] + self.X['relax_time'] + 1e-6)
        
        # 3. 咖啡因摄入率 - 衡量每工作小时摄入的咖啡因量
        self.X['caffeine_per_hour'] = self.X['caffeine_intake'] / (self.X['work_time'] + 1e-6)
        
        # 4. 屏幕时间比例 - 评估手机使用在日常活动中的占比
        self.X['screen_time_ratio'] = self.X['phone_time'] / (self.X['total_activities'] + 1e-6)
        
        # 5. 健康活动指数 - 较高的值表示更健康的活动模式
        self.X['healthy_activity_index'] = (self.X['exercise_time'] + self.X['reading_time']) / (self.X['phone_time'] + 1e-6)
        
        # 6. 每日总时间使用 - 评估一天中时间分配的总和
        self.X['daily_time_utilization'] = (self.X['work_time'] + self.X['total_activities'])
        
        print(f"特征工程后的数据集形状: {self.X.shape}")
        print("新创建的特征:")
        new_features = ['total_activities', 'work_life_balance', 'caffeine_per_hour', 
                       'screen_time_ratio', 'healthy_activity_index', 'daily_time_utilization']
        print(self.X[new_features].describe())
        
        # 绘制新特征与目标变量的相关性
        print("\n新特征与睡眠时间的相关性:")
        new_corr = self.X[new_features].corrwith(self.y).sort_values(ascending=False)
        for feat, corr in new_corr.items():
            print(f"{feat}: {corr:.3f}")
        
        return self
    
    def prepare_data(self, test_size=0.2, random_state=42):
        """
        准备训练和测试数据
        
        此方法执行两项关键任务:
        1. 将数据集随机分割为训练集和测试集
        2. 对特征进行标准化处理，使不同量级的特征具有相同的尺度
        
        标准化可以显著提高某些算法（如线性模型、SVM等）的性能，
        确保特征值在相似范围内，避免某些特征因数值较大而主导模型训练。
        
        参数:
            test_size (float): 测试集比例，默认为0.2（即20%的数据用于测试）
            random_state (int): 随机种子，确保结果可复现
            
        返回:
            self: 返回实例本身以支持方法链式调用
        """
        print("\n=== 数据准备 ===")
        
        # 分割数据集 - 使用分层抽样确保训练集和测试集的分布相似
        self.X_train, self.X_test, self.y_train, self.y_test = train_test_split(
            self.X, self.y, test_size=test_size, random_state=random_state
        )
        
        # 标准化特征 - 在训练集上拟合标准化器并转换
        # 对测试集应用相同的变换，确保数据处理的一致性
        self.X_train_scaled = self.scaler.fit_transform(self.X_train)
        self.X_test_scaled = self.scaler.transform(self.X_test)
        
        print(f"训练集大小: {self.X_train.shape[0]}")
        print(f"测试集大小: {self.X_test.shape[0]}")
        
        return self
    
    def train_models(self):
        """
        训练多个回归模型并评估性能
        
        此方法实现了以下功能:
        1. 初始化多种不同类型的回归算法模型
        2. 针对不同类型的模型使用不同的数据预处理方式
        3. 训练每个模型并评估其性能指标
        4. 确定最佳性能模型并进行可视化比较
        
        模型选择策略:
        - 线性模型(线性回归、岭回归、Lasso回归): 使用标准化后的数据，有助于提高收敛速度和精度
        - 树模型(随机森林、梯度提升): 使用原始特征数据，不受特征尺度影响
        
        评估指标包括均方误差(MSE)、均方根误差(RMSE)、平均绝对误差(MAE)、
        R²分数、调整R²、平均绝对百分比误差(MAPE)和对称平均绝对百分比误差(SMAPE)。
        
        返回:
            self: 返回实例本身以支持方法链式调用
        """
        print("\n=== 模型训练 ===")
        
        # 定义模型
        models = {
            '线性回归': LinearRegression(),
            '岭回归': Ridge(random_state=42),
            'Lasso回归': Lasso(random_state=42),
            '随机森林': RandomForestRegressor(random_state=42),
            '梯度提升': GradientBoostingRegressor(random_state=42)
        }
        
        # 训练和评估模型
        for name, model in models.items():
            # 根据模型类型选择数据预处理方式
            if name in ['线性回归', '岭回归', 'Lasso回归']:
                # 线性模型需要标准化特征
                model.fit(self.X_train_scaled, self.y_train)
                y_pred = model.predict(self.X_test_scaled)
            else:
                # 树模型不需要标准化
                model.fit(self.X_train, self.y_train)
                y_pred = model.predict(self.X_test)
            
            # 评估模型并保存结果
            self.results[name] = self._evaluate_model(self.y_test, y_pred, name)
        
        # 找出最佳模型（基于R²分数）
        best_model_name = max(self.results, key=lambda x: self.results[x]['R2'])
        self.best_model_name = best_model_name
        
        print(f"\n最佳模型: {best_model_name} (R² = {self.results[best_model_name]['R2']:.4f})")
        
        # 可视化模型性能比较
        self._plot_model_comparison()
        
        return self
    
    def _evaluate_model(self, y_true, y_pred, model_name):
        """
        评估模型性能并返回详细指标
        
        计算并返回模型的全面性能指标，包括基础误差指标和高级统计指标
        这些指标从不同角度评估模型预测能力，帮助全面了解模型表现。
        
        参数:
            y_true (array-like): 真实目标值（实际睡眠时间）
            y_pred (array-like): 模型预测值
            model_name (str): 模型名称，用于输出标识
            
        返回:
            dict: 包含所有评估指标的字典
                - MSE: 均方误差
                - RMSE: 均方根误差（小时）
                - MAE: 平均绝对误差（小时）
                - R2: 决定系数
                - 调整R²: 考虑特征数量的决定系数
                - MAPE: 平均绝对百分比误差（%）
                - SMAPE: 对称平均绝对百分比误差（%）
        """
        # 基础评估指标
        mse = mean_squared_error(y_true, y_pred)  # 均方误差 - 平方放大了大误差的影响
        rmse = np.sqrt(mse)  # 均方根误差 - 与目标变量单位一致，更易解释
        mae = mean_absolute_error(y_true, y_pred)  # 平均绝对误差 - 直观的误差度量
        r2 = r2_score(y_true, y_pred)  # 决定系数 - 模型解释目标变量变异的能力
        
        # 计算调整后的R² - 考虑了特征数量的影响，更适合比较不同特征集的模型
        n = len(y_true)  # 样本数量
        p = self.X_train.shape[1]  # 特征数量
        adjusted_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)
        
        # 计算MAPE和SMAPE - 相对误差指标，更易理解百分比误差
        mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100  # 平均绝对百分比误差
        smape = np.mean(2 * np.abs(y_pred - y_true) / (np.abs(y_pred) + np.abs(y_true))) * 100  # 对称平均绝对百分比误差
        
        # 输出详细评估结果
        print(f"\n{model_name} 评估结果:")
        print(f"均方误差 (MSE): {mse:.4f}")
        print(f"均方根误差 (RMSE): {rmse:.4f} 小时")
        print(f"平均绝对误差 (MAE): {mae:.4f} 小时")
        print(f"R² 分数: {r2:.4f}")
        print(f"调整R²: {adjusted_r2:.4f}")
        print(f"平均绝对百分比误差 (MAPE): {mape:.2f}%")
        print(f"对称平均绝对百分比误差 (SMAPE): {smape:.2f}%")
        
        # 返回所有评估指标
        return {
            'MSE': mse, 
            'RMSE': rmse, 
            'MAE': mae, 
            'R2': r2,
            '调整R²': adjusted_r2,
            'MAPE': mape,
            'SMAPE': smape
        }
    
    def _plot_learning_curve(self):
        """
        绘制学习曲线，展示模型在不同训练集大小下的性能
        """
        try:
            from sklearn.model_selection import learning_curve
            
            # 定义训练集大小范围
            train_sizes = np.linspace(0.1, 1.0, 10)
            
            # 根据模型类型选择数据
            if self.best_model_name in ['线性回归', '岭回归', 'Lasso回归']:
                X = self.X_train_scaled
                y = self.y_train
            else:
                X = self.X_train
                y = self.y_train
            
            # 计算学习曲线
            train_sizes, train_scores, validation_scores = learning_curve(
                estimator=self.best_model,
                X=X,
                y=y,
                train_sizes=train_sizes,
                cv=5,
                scoring='r2',
                n_jobs=-1
            )
            
            # 计算均值和标准差
            train_mean = np.mean(train_scores, axis=1)
            train_std = np.std(train_scores, axis=1)
            validation_mean = np.mean(validation_scores, axis=1)
            validation_std = np.std(validation_scores, axis=1)
            
            # 绘制学习曲线
            plt.figure(figsize=(12, 8))
            plt.plot(train_sizes, train_mean, color='blue', marker='o', 
                    markersize=5, label='训练集 R²')
            plt.fill_between(train_sizes, train_mean + train_std, train_mean - train_std,
                            alpha=0.15, color='blue')
            plt.plot(train_sizes, validation_mean, color='green', marker='s', 
                    markersize=5, label='验证集 R²')
            plt.fill_between(train_sizes, validation_mean + validation_std, 
                            validation_mean - validation_std, alpha=0.15, color='green')
            
            # 分析欠拟合或过拟合
            train_final = train_mean[-1]
            val_final = validation_mean[-1]
            gap = abs(train_final - val_final)
            
            # 添加诊断信息
            if train_final < 0.5 and val_final < 0.5:
                plt.axvline(x=train_sizes[-1], color='red', linestyle='--', alpha=0.7)
                plt.text(train_sizes[-1] + 0.01, 0.5, '欠拟合', color='red')
            elif gap > 0.2:
                plt.axvline(x=train_sizes[-1], color='orange', linestyle='--', alpha=0.7)
                plt.text(train_sizes[-1] + 0.01, 0.5, '过拟合', color='orange')
            else:
                plt.axvline(x=train_sizes[-1], color='green', linestyle='--', alpha=0.7)
                plt.text(train_sizes[-1] + 0.01, 0.5, '良好拟合', color='green')
            
            plt.title(f'{self.best_model_name} 学习曲线')
            plt.xlabel('训练样本数')
            plt.ylabel('R² 分数')
            plt.legend(loc='lower right')
            plt.grid(True, alpha=0.3)
            plt.ylim(-0.1, 1.1)
            plt.tight_layout()
            plt.savefig('学习曲线.png', dpi=300, bbox_inches='tight')
            plt.show()
            
            # 输出诊断结果
            print(f"学习曲线分析: 训练集R²={train_final:.4f}, 验证集R²={val_final:.4f}, 差距={gap:.4f}")
            if train_final < 0.5 and val_final < 0.5:
                print("诊断: 模型可能存在欠拟合。建议增加模型复杂度或添加更多特征。")
            elif gap > 0.2:
                print("诊断: 模型可能存在过拟合。建议增加正则化或使用更多训练数据。")
            else:
                print("诊断: 模型拟合状态良好。")
                
        except Exception as e:
            print(f"无法绘制学习曲线: {e}")
    
    def _plot_prediction_distribution(self):
        """
        绘制预测值分布和真实值分布的对比直方图
        """
        try:
            # 获取预测值
            if self.best_model_name in ['线性回归', '岭回归', 'Lasso回归']:
                y_pred = self.best_model.predict(self.X_test_scaled)
            else:
                y_pred = self.best_model.predict(self.X_test)
            
            # 绘制对比直方图
            plt.figure(figsize=(12, 6))
            bins = 20
            
            plt.hist(self.y_test, bins=bins, alpha=0.5, color='blue', label='真实睡眠时间')
            plt.hist(y_pred, bins=bins, alpha=0.5, color='green', label='预测睡眠时间')
            
            plt.xlabel('睡眠时间 (小时)')
            plt.ylabel('频率')
            plt.title(f'{self.best_model_name} 真实值与预测值分布对比')
            plt.legend(loc='upper right')
            plt.grid(True, alpha=0.3)
            plt.tight_layout()
            plt.savefig('预测分布对比.png', dpi=300, bbox_inches='tight')
            plt.show()
            
        except Exception as e:
            print(f"无法绘制预测分布: {e}")
    
    def _plot_model_comparison(self):
        """可视化模型比较"""
        results_df = pd.DataFrame(self.results).T
        
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        
        # RMSE 比较
        axes[0,0].bar(results_df.index, results_df['RMSE'], color='skyblue')
        axes[0,0].set_title('RMSE 比较')
        axes[0,0].set_ylabel('RMSE')
        for i, v in enumerate(results_df['RMSE']):
            axes[0,0].text(i, v, f'{v:.3f}', ha='center', va='bottom')
        
        # MAE 比较
        axes[0,1].bar(results_df.index, results_df['MAE'], color='lightgreen')
        axes[0,1].set_title('MAE 比较')
        axes[0,1].set_ylabel('MAE')
        for i, v in enumerate(results_df['MAE']):
            axes[0,1].text(i, v, f'{v:.3f}', ha='center', va='bottom')
        
        # R² 比较
        axes[1,0].bar(results_df.index, results_df['R2'], color='lightcoral')
        axes[1,0].set_title('R² 比较')
        axes[1,0].set_ylabel('R²')
        for i, v in enumerate(results_df['R2']):
            axes[1,0].text(i, v, f'{v:.3f}', ha='center', va='bottom')
        
        # MSE 比较
        axes[1,1].bar(results_df.index, results_df['MSE'], color='lightyellow')
        axes[1,1].set_title('MSE 比较')
        axes[1,1].set_ylabel('MSE')
        for i, v in enumerate(results_df['MSE']):
            axes[1,1].text(i, v, f'{v:.3f}', ha='center', va='bottom')
        
        plt.tight_layout()
        plt.savefig('模型性能比较.png', dpi=300, bbox_inches='tight')
        plt.show()
        
        # 雷达图比较所有指标
        if 'MAPE' in results_df.columns:
            metrics = ['RMSE', 'MAE', 'MSE', 'R2', 'MAPE']
            # 标准化指标（除了R2）
            radar_df = results_df.copy()
            for metric in ['RMSE', 'MAE', 'MSE', 'MAPE']:
                radar_df[metric] = 1 - radar_df[metric] / radar_df[metric].max()
            
            angles = np.linspace(0, 2*np.pi, len(metrics), endpoint=False).tolist()
            angles += angles[:1]  # 闭合雷达图
            
            plt.figure(figsize=(10, 8))
            ax = plt.subplot(111, polar=True)
            
            for model in radar_df.index:
                values = radar_df.loc[model, metrics].tolist()
                values += values[:1]  # 闭合雷达图
                ax.plot(angles, values, linewidth=2, label=model)
                ax.fill(angles, values, alpha=0.1)
            
            ax.set_thetagrids(np.degrees(angles[:-1]), metrics)
            ax.set_ylim(0, 1)
            plt.legend(loc='upper right', bbox_to_anchor=(0.1, 0.1))
            plt.title('模型性能雷达图比较')
            plt.tight_layout()
            plt.savefig('模型性能雷达图.png', dpi=300, bbox_inches='tight')
            plt.show()
        
        print("\n模型性能总结:")
        print(results_df)
    
    def optimize_best_model(self):
        """
        优化最佳模型的超参数
        """
        print("\n=== 模型超参数优化 ===")
        
        # 根据最佳模型类型进行超参数调优
        if self.best_model_name == '随机森林':
            param_grid = {
                'n_estimators': [50, 100, 200],
                'max_depth': [None, 10, 20, 30],
                'min_samples_split': [2, 5, 10],
                'min_samples_leaf': [1, 2, 4]
            }
            self.best_model = RandomForestRegressor(random_state=42)
            
        elif self.best_model_name == '梯度提升':
            param_grid = {
                'n_estimators': [50, 100, 200],
                'learning_rate': [0.01, 0.1, 0.2],
                'max_depth': [3, 5, 7],
                'subsample': [0.8, 1.0]
            }
            self.best_model = GradientBoostingRegressor(random_state=42)
            
        elif self.best_model_name == '岭回归':
            param_grid = {
                'alpha': [0.01, 0.1, 1.0, 10.0, 100.0]
            }
            self.best_model = Ridge(random_state=42)
            
        elif self.best_model_name == 'Lasso回归':
            param_grid = {
                'alpha': [0.01, 0.1, 1.0, 10.0, 100.0]
            }
            self.best_model = Lasso(random_state=42)
            
        else:  # 线性回归
            print("线性回归无需超参数优化")
            self.best_model = LinearRegression()
            self.best_model.fit(self.X_train_scaled, self.y_train)
            return self
        
        # 使用网格搜索进行超参数调优
        grid_search = GridSearchCV(estimator=self.best_model, param_grid=param_grid,
                                 cv=5, scoring='r2', n_jobs=-1, verbose=2)
        
        if self.best_model_name in ['线性回归', '岭回归', 'Lasso回归']:
            grid_search.fit(self.X_train_scaled, self.y_train)
        else:
            grid_search.fit(self.X_train, self.y_train)
        
        # 获取最佳参数
        best_params = grid_search.best_params_
        self.best_model = grid_search.best_estimator_
        
        print(f"\n{self.best_model_name} 最佳参数:")
        for param, value in best_params.items():
            print(f"  {param}: {value}")
        
        # 使用最佳模型进行预测
        if self.best_model_name in ['线性回归', '岭回归', 'Lasso回归']:
            y_pred_optimized = self.best_model.predict(self.X_test_scaled)
        else:
            y_pred_optimized = self.best_model.predict(self.X_test)
        
        # 评估优化后的模型
        optimized_results = self._evaluate_model(self.y_test, y_pred_optimized, 
                                               f"{self.best_model_name} (优化后)")
        
        # 交叉验证
        print("\n=== 交叉验证 ===")
        if self.best_model_name in ['线性回归', '岭回归', 'Lasso回归']:
            cv_scores = cross_val_score(self.best_model, self.X_train_scaled, self.y_train, 
                                      cv=5, scoring='r2')
        else:
            cv_scores = cross_val_score(self.best_model, self.X_train, self.y_train, 
                                      cv=5, scoring='r2')
        
        print(f"5折交叉验证 R² 分数: {cv_scores.mean():.4f} ± {cv_scores.std():.4f}")
        
        return self
    
    def analyze_results(self):
        """
        分析模型结果和特征重要性
        """
        print("\n=== 模型结果分析 ===")
        
        # 预测结果可视化
        if self.best_model_name in ['线性回归', '岭回归', 'Lasso回归']:
            y_pred_best = self.best_model.predict(self.X_test_scaled)
        else:
            y_pred_best = self.best_model.predict(self.X_test)
        
        # 基础预测结果可视化
        self._plot_prediction_results(y_pred_best)
        
        # 高级可视化和分析
        print("\n=== 高级模型诊断 ===")
        
        # 学习曲线分析
        self._plot_learning_curve()
        
        # 预测分布对比
        self._plot_prediction_distribution()
        
        # 特征重要性分析（如果模型支持）
        self._analyze_feature_importance()
        
        # 对于随机森林模型绘制验证曲线
        if '随机森林' in self.best_model_name:
            self._plot_validation_curve()
        
        return self
    
    def _plot_prediction_results(self, y_pred):
        """绘制预测结果可视化"""
        plt.figure(figsize=(14, 6))
        
        # 预测值 vs 实际值
        plt.subplot(1, 2, 1)
        plt.scatter(self.y_test, y_pred, alpha=0.6, color='blue')
        plt.plot([self.y_test.min(), self.y_test.max()], [self.y_test.min(), self.y_test.max()], 'r--', lw=2)
        plt.xlabel('实际睡眠时间')
        plt.ylabel('预测睡眠时间')
        plt.title(f'{self.best_model_name}：预测值 vs 实际值')
        plt.grid(True, alpha=0.3)
        
        # 残差图
        plt.subplot(1, 2, 2)
        residuals = self.y_test - y_pred
        plt.scatter(y_pred, residuals, alpha=0.6, color='green')
        plt.axhline(y=0, color='r', linestyle='--')
        plt.xlabel('预测睡眠时间')
        plt.ylabel('残差')
        plt.title(f'{self.best_model_name}：残差图')
        plt.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.savefig('预测结果可视化.png', dpi=300, bbox_inches='tight')
        plt.show()
        
        # 残差分布
        plt.figure(figsize=(8, 6))
        plt.hist(residuals, bins=30, alpha=0.7, color='purple')
        plt.xlabel('残差')
        plt.ylabel('频次')
        plt.title(f'{self.best_model_name}：残差分布')
        plt.grid(True, alpha=0.3)
        plt.savefig('残差分布.png', dpi=300, bbox_inches='tight')
        plt.show()
    
    def _analyze_feature_importance(self):
        """分析特征重要性"""
        try:
            # 对于树模型
            if hasattr(self.best_model, 'feature_importances_'):
                importances = self.best_model.feature_importances_
                indices = np.argsort(importances)[::-1]
                
                # 创建水平条形图
                plt.figure(figsize=(14, 10))
                bars = plt.barh(range(len(indices)), importances[indices], color='lightblue')
                plt.yticks(range(len(indices)), self.X.columns[indices], fontsize=12)
                plt.xlabel('重要性', fontsize=14)
                plt.title(f'{self.best_model_name} 特征重要性', fontsize=16)
                
                # 添加重要性数值标签
                for i, bar in enumerate(bars):
                    plt.text(bar.get_width() + 0.005, bar.get_y() + bar.get_height()/2, 
                           f'{importances[indices[i]]:.4f}', va='center', fontsize=10)
                
                plt.tight_layout()
                plt.savefig('特征重要性.png', dpi=300, bbox_inches='tight')
                plt.show()
                
                # 特征重要性累积图
                plt.figure(figsize=(12, 6))
                cumulative_importance = np.cumsum(sorted(importances, reverse=True))
                plt.plot(range(1, len(cumulative_importance) + 1), cumulative_importance, 
                        'o-', color='red')
                plt.axhline(y=0.95, color='gray', linestyle='--', alpha=0.7)
                plt.xlabel('特征数量')
                plt.ylabel('累积重要性')
                plt.title('特征重要性累积曲线')
                plt.grid(True, alpha=0.3)
                plt.tight_layout()
                plt.savefig('特征重要性累积曲线.png', dpi=300, bbox_inches='tight')
                plt.show()
                
                print("\n特征重要性排序:")
                for i in indices:
                    print(f"{self.X.columns[i]}: {importances[i]:.4f}")
                    
                # 计算达到95%重要性需要的特征数量
                n_features_95 = np.argmax(cumulative_importance >= 0.95) + 1
                print(f"\n达到95%重要性需要前{n_features_95}个特征")
            
            # 对于线性模型
            elif hasattr(self.best_model, 'coef_'):
                coefficients = self.best_model.coef_
                
                # 按绝对值排序
                indices = np.argsort(np.abs(coefficients))[::-1]
                
                # 创建水平条形图
                plt.figure(figsize=(14, 10))
                # 使用不同颜色表示正负系数
                colors = ['red' if coefficients[i] < 0 else 'blue' for i in indices]
                bars = plt.barh(range(len(indices)), np.abs(coefficients)[indices], color=colors)
                plt.yticks(range(len(indices)), self.X.columns[indices], fontsize=12)
                plt.xlabel('系数绝对值', fontsize=14)
                plt.title(f'{self.best_model_name} 特征系数（按绝对值排序）', fontsize=16)
                
                # 添加系数数值标签
                for i, bar in enumerate(bars):
                    actual_coef = coefficients[indices[i]]
                    plt.text(bar.get_width() + 0.005, bar.get_y() + bar.get_height()/2, 
                           f'{actual_coef:.4f}', va='center', fontsize=10)
                
                # 添加图例
                plt.legend(['负相关', '正相关'], loc='lower right')
                
                plt.tight_layout()
                plt.savefig('特征系数.png', dpi=300, bbox_inches='tight')
                plt.show()
                
                # 特征系数分布图
                plt.figure(figsize=(12, 6))
                plt.bar(range(len(coefficients)), coefficients, color=colors)
                plt.axhline(y=0, color='black', linestyle='-', linewidth=0.8)
                plt.xticks(range(len(coefficients)), self.X.columns, rotation=90, fontsize=10)
                plt.xlabel('特征')
                plt.ylabel('系数值')
                plt.title('各特征系数分布')
                plt.grid(True, alpha=0.3, axis='y')
                plt.tight_layout()
                plt.savefig('特征系数分布.png', dpi=300, bbox_inches='tight')
                plt.show()
                
                print("\n特征系数排序:")
                for i in indices:
                    print(f"{self.X.columns[i]}: {coefficients[i]:.4f}")
        
        except Exception as e:
            print(f"无法计算特征重要性: {e}")
    
    def _plot_validation_curve(self):
        """
        绘制验证曲线（针对n_estimators参数）
        """
        try:
            from sklearn.model_selection import validation_curve
            
            # 定义参数范围
            param_range = [10, 50, 100, 150, 200, 250, 300]
            
            # 根据模型类型选择数据
            if self.best_model_name in ['线性回归', '岭回归', 'Lasso回归']:
                X = self.X_train_scaled
                y = self.y_train
            else:
                X = self.X_train
                y = self.y_train
            
            # 计算验证曲线
            train_scores, test_scores = validation_curve(
                self.best_model, X, y, param_name='n_estimators', 
                param_range=param_range, cv=5, scoring='r2', n_jobs=-1
            )
            
            # 计算均值和标准差
            train_mean = np.mean(train_scores, axis=1)
            train_std = np.std(train_scores, axis=1)
            test_mean = np.mean(test_scores, axis=1)
            test_std = np.std(test_scores, axis=1)
            
            # 绘制验证曲线
            plt.figure(figsize=(12, 8))
            plt.plot(param_range, train_mean, color='blue', marker='o', 
                    markersize=5, label='训练集 R²')
            plt.fill_between(param_range, train_mean + train_std, train_mean - train_std,
                            alpha=0.15, color='blue')
            plt.plot(param_range, test_mean, color='green', marker='s', 
                    markersize=5, label='验证集 R²')
            plt.fill_between(param_range, test_mean + test_std, test_mean - test_std,
                            alpha=0.15, color='green')
            
            # 标记最佳参数
            best_n_estimators = param_range[np.argmax(test_mean)]
            best_score = np.max(test_mean)
            plt.scatter(best_n_estimators, best_score, color='red', s=200, 
                       zorder=5, label=f'最佳参数: {best_n_estimators}')
            
            plt.title(f'{self.best_model_name} 验证曲线 (n_estimators)')
            plt.xlabel('n_estimators')
            plt.ylabel('R² 分数')
            plt.legend(loc='lower right')
            plt.grid(True, alpha=0.3)
            plt.ylim(-0.1, 1.1)
            plt.savefig('验证曲线.png', dpi=300, bbox_inches='tight')
            plt.show()
            
            print(f"验证曲线分析: 最佳n_estimators参数为{best_n_estimators}，对应的验证R²为{best_score:.4f}")
            
        except Exception as e:
            print(f"无法绘制验证曲线: {e}")
    
    def generate_conclusions(self):
        """
        生成结论和建议
        """
        print("\n=== 结论和建议 ===")
        
        print("\n1. 数据概况:")
        print(f"   - 原始数据集包含 {self.original_shape[0]} 条记录，{self.original_shape[1]} 个特征")
        print(f"   - 处理后数据集包含 {len(self.df)} 条记录")
        print(f"   - 最终特征数: {len(self.X.columns)} 个")
        print(f"   - 睡眠时间范围: {self.df['sleep_time'].min():.1f} - {self.df['sleep_time'].max():.1f} 小时")
        print(f"   - 平均睡眠时间: {self.df['sleep_time'].mean():.1f} 小时")
        
        print("\n2. 关键发现:")
        # 根据相关性分析生成关键发现
        correlation_with_sleep = self.df.corr()['sleep_time'].drop('sleep_time')
        
        # 正相关特征（对睡眠有帮助）
        positive_corr = correlation_with_sleep[correlation_with_sleep > 0].sort_values(ascending=False)
        if not positive_corr.empty:
            top_positive = positive_corr.index[0] if len(positive_corr) > 0 else None
            if top_positive:
                print(f"   - {top_positive.replace('_', ' ')}与睡眠时间呈正相关")
        
        # 负相关特征（对睡眠有负面影响）
        negative_corr = correlation_with_sleep[correlation_with_sleep < 0].sort_values()
        if not negative_corr.empty:
            top_negative = negative_corr.index[0] if len(negative_corr) > 0 else None
            if top_negative:
                print(f"   - {top_negative.replace('_', ' ')}与睡眠时间呈负相关")
        
        # 通用发现
        print("   - 运动时间通常对睡眠时间有积极影响")
        print("   - 工作时间和咖啡因摄入量通常对睡眠时间有负面影响")
        print("   - 放松时间有助于增加睡眠时间")
        
        best_model_results = self.results.get(self.best_model_name, {})
        print("\n3. 模型性能:")
        print(f"   - 最佳模型: {self.best_model_name}")
        print(f"   - R² 分数: {best_model_results.get('R2', 0):.3f}")
        print(f"   - 预测误差: RMSE = {best_model_results.get('RMSE', 0):.3f} 小时")
        print(f"   - 平均绝对误差: MAE = {best_model_results.get('MAE', 0):.3f} 小时")
        
        print("\n4. 健康建议:")
        print("   - 增加日常运动量，尤其是有氧运动")
        print("   - 减少咖啡因摄入，特别是在下午和晚上")
        print("   - 控制工作时间，保持良好的工作生活平衡")
        print("   - 减少睡前手机等电子设备的使用时间")
        print("   - 确保充足的放松和休闲时间")
        print("   - 建立规律的睡眠习惯")
        
        print("\n5. 模型应用:")
        print("   - 可用于健康应用程序的睡眠预测功能")
        print("   - 为个人提供个性化的睡眠改善建议")
        print("   - 帮助识别影响睡眠的关键生活方式因素")
        print("   - 作为睡眠研究的辅助工具")
        
        return self
    
    def save_model(self, model_path='sleep_prediction_model.pkl', scaler_path='sleep_scaler.pkl'):
        """
        保存训练好的模型和标准化器
        
        参数:
            model_path: 模型保存路径
            scaler_path: 标准化器保存路径
        """
        # 保存模型和标准化器
        joblib.dump(self.best_model, model_path)
        joblib.dump(self.scaler, scaler_path)
        
        print(f"\n模型已保存为 '{model_path}'")
        print(f"标准化器已保存为 '{scaler_path}'")
        
        # 创建使用示例
        self._create_usage_example()
        
        return self
    
    def _create_usage_example(self):
        """创建模型使用示例"""
        print("\n=== 使用示例 ===")
        print("# 加载模型和标准化器")
        print("import joblib")
        print("import numpy as np")
        print("import pandas as pd")
        print()
        print("model = joblib.load('sleep_prediction_model.pkl')")
        print("scaler = joblib.load('sleep_scaler.pkl')")
        print()
        print("# 创建预测函数")
        print("def predict_sleep_time(exercise_time, reading_time, phone_time, work_time, caffeine_intake, relax_time):")
        print("    # 创建输入数据")
        print("    input_data = pd.DataFrame({")
        print("        'exercise_time': [exercise_time],")
        print("        'reading_time': [reading_time],")
        print("        'phone_time': [phone_time],")
        print("        'work_time': [work_time],")
        print("        'caffeine_intake': [caffeine_intake],")
        print("        'relax_time': [relax_time]")
        print("    })")
        print()
        print("    # 创建与训练时相同的衍生特征")
        print("    input_data['total_activities'] = input_data['exercise_time'] + input_data['reading_time'] + ")
        print("                                  input_data['phone_time'] + input_data['relax_time']")
        print("    input_data['work_life_balance'] = input_data['work_time'] / (input_data['exercise_time'] + input_data['relax_time'] + 1e-6)")
        print("    input_data['caffeine_per_hour'] = input_data['caffeine_intake'] / (input_data['work_time'] + 1e-6)")
        print("    input_data['screen_time_ratio'] = input_data['phone_time'] / (input_data['total_activities'] + 1e-6)")
        print("    input_data['healthy_activity_index'] = (input_data['exercise_time'] + input_data['reading_time']) / (input_data['phone_time'] + 1e-6)")
        print("    input_data['daily_time_utilization'] = input_data['work_time'] + input_data['total_activities']")
        print()
        print("    # 根据模型类型选择是否需要标准化")
        print("    if model.__class__.__name__ in ['LinearRegression', 'Ridge', 'Lasso']:")
        print("        scaled_data = scaler.transform(input_data)")
        print("        prediction = model.predict(scaled_data)")
        print("    else:")
        print("        prediction = model.predict(input_data)")
        print()
        print("    return prediction[0]")
        print()
        print("# 示例使用")
        print("sleep_hours = predict_sleep_time(1.5, 0.8, 2.0, 8.0, 150, 1.2)")
        print("print(f'预测睡眠时间: {sleep_hours:.2f} 小时')")

# 主函数
def main():
    print("=== 睡眠时间预测机器学习项目 ===\n")
    
    # 初始化系统
    sleep_system = SleepPredictionSystem()
    
    try:
        # 执行完整的工作流程
        sleep_system.load_data()\
                   .handle_outliers(method='iqr')\
                   .exploratory_data_analysis()\
                   .feature_engineering()\
                   .prepare_data()\
                   .train_models()\
                   .optimize_best_model()\
                   .analyze_results()\
                   .generate_conclusions()\
                   .save_model()
        
        print("\n=== 项目完成 ===")
        print("所有步骤已成功执行！")
        
    except Exception as e:
        print(f"\n发生错误: {e}")
        import traceback
        traceback.print_exc()

if __name__ == "__main__":
    main()